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ABSTRACT 

In this paper we investigate the space and velocity distributions of old neutron stars (aged 
10 9 to 10 10 yr) in our Galaxy. Galactic old Neutron Stars (NSs) population fills a torus-like area 
extending to a few tens kiloparsccs above the galactic plane. The initial velocity distribution of 
NSs is not well known, in this work we adopt a three component initial distribution, as given by 
the contribution of kick velocities, circular velocities and Maxwellian velocities. For the spatial 
initial distribution we use a T function. We then use Monte Carlo simulations to follow the 
evolution of the NSs under the influence of the Paczyhski Galactic gravitational potential. Our 
calculations show that NS orbits have a very large Galactic radial expansion and that their 
radial distribution peak is quite close to their progenitors' one. We also study the NS vertical 
distribution and find that it can well be described by a double exponential low. Finally, we 
investigate the correlation of the vertical and radial distribution and study the radial dependence 
of scale-heights. 
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Introduction 
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It is commonly accepted that Neutron Stars 
(NSs) are born when massive OB-stars exhaust 
their nuclear fuel and end their lives in core- 
collapse supern ova explosions, near the Galactic 
disk ( see e.g. iBhattacharva fc van den Heuvel 



and that they have then moved away from 
the Galactic plane with ave rage kick velocities 
around 200 - 500 km/s (e. g. ICordes fc Chernoff I 
(<1998h . iLvne et all (<2004h . iHobbs et all (12005^ . 
About 10 9 NSs are thought to populate our 
Galaxy, but only 2 x 10 3 are directly observed as 
radi o pulsars or as accr etion-powered Xray bina- 
ries ( Sartore et al.l [2010h . as a consequence, little 
is known about their statistical properties. The 
estimation of pulsar velocities relies on direct dis- 
tance measurements, which can be obtained by 
dispersion measures tog ether with a Galactic elec - 
tron density model (e. g. iTavlor fc Cordes ( 19931 ). 
Cordes fc La"ziol(|2002l) V The mechanis ms for pro 



ducing high velocities are still unknown ([Lorimer 
2008h . 



Numerical simulations are valuable tools for un- 
derstanding the spatial and velocity distribution 
of NSs in the Galaxy, and they have been used 
by several authors ( see e.g. ICaldwell fc Ostriker 
(<198ll ) . lOfek I (l2009h . iKatsanikas fc Patsis I (|20 1 lh ) . 
In particular, Paczyhski ( 1990h (hereafter P90) 
simulated the motion of NSs in connection with 
the galactic origin of gamma ray bursts and cal- 
culated the NS space density distribution. In 
the same work Paczyhski also suggested a sim- 
plified expression for the gravitational potential 
which is still often applied in the simul a tion o f 
NS distribution. For example IWei et al.1 (|2010ah 
used this potential. They considered the old 
NSs, i.e. NSs whose age is between 10 9 and 
10 10 yr. They adopted a Galactic distribution 
with one-component initial random velocity mod- 
els. The aim of the pres ent work i s to imp rove 
the model developed in I Wei et al.l (|2010al ) for 
studying the distribution of old NSs as a func- 
tion of the initial position distribution, of the 
initial velocity distribution and of the galactic 
gravitational potential. Based on P90 gravita- 
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tional potential we consider the evolution of a two- 
component Maxwellian initial random velocity dis- 
tribution (we ado pt the velocity distr i butio ns of 
NSs at birth from Arzoumanian et al. I d2002h and 



Faucher-Giguere fe Kaspi ( 20061 )). We perform 



integration of NS velocities using Monte Carlo in- 
tegration techniques with different conditions de- 
veloped for this purpose. We also aim at obtaining 
the NS trajectories under variety of assumptions. 

The paper is structured as follows: in Sec. [2] we 
describe the ingredients of the simulation, i.e. the 
NS initial position and velocity distribution and 
the Galactic gravitational potential. We present 
the results of the simulation in Sec. [3] and |U In 
Sec.[3]we investigate the NSs orbits; while in Sec. 2] 
we investigate their vertical and radial distribu- 
tions. We fit the two exponential decay model, for 
each spacing segment of R to derive the high scale 
heights. We fit also the R distributions at differ- 
ent scale heights. We discuss our results and their 
possible implications in Sec. [5j 

2. Simulation ingredients 

In this Section we present the ingredients of the 
model: the NS birth rate and Monte Carlo sim- 
ulation; the gravitational potential; the position 
distributions of NSs; and the initial velocities and 
equations of motion. 

2.1. NSs birth rate and Monte Carlo sim- 
ulation 

Theoretically, estimating the birthrate of a pop- 
ulation of sources is simple. However, for the NS 
population, precise estimates of both the num- 
ber and lifetime of the sources are hard to ob- 
tain, because they may have been heavily bi- 
ased by a number of observational selection ef- 
fects. The birthrate of NSs (77) within the whole 
Galactic disk is roughly 0.9 ~ 2 per cent ury 
(|Lvne fe Graham-Smith ll2007tlLorimer 1120081 ). If 



one assumes a life time r ~ 10 9 - 10 10 yr, then an 
estimate of the total number of NSs is: 



N ■■ 



t * r\ 



10 



7-8 



(1) 



The problem is best tackled using a Monte 
Carlo simulations of NS positions, orbits and ve- 
locities, taking into account the birthplaces and 



the initial kick velocities. We study the resulting 
phase spatial distributions concerning the Galactic 
potential and the distribution of progenitors and 
birth velocities, focusing on the numerical prop- 
erties of the NS populations in the disk and in 
the solar neighborhood. NS orbits are obtained 
by solving the equations of motion in the P90's 
gravitational potential. 

2.2. Galactic gravitational potential 

It is known that the Galactic gravitational po- 
tential causes oscillation of objects along the di- 
recti on perpendicular t o the Galactic plane (see 



e.g. iLvne et al. I (|1982l ). and references therein). 
To track the evolution and motion of NSs popula- 
tion, the gravitational potential P90 is taken to be 
a homogeneous function of the density, and ignore 
the interstellar friction. This is a reliable approx- 
imation for our axisymmetric model, because the 
steady state distribution of old NSs depends only 
weakly on th e non- homogeneo us part of the galac- 
tic potential ( Frei et al1ll992l ). Asymmetry in the 
kinematics, which is likely due to the finite lifetime 
of the stars and Galac tic potential struct ure, is a 



20091 ). How- 



rclatively small effect (jPerets et al. 
ever, using P90 may not be a good approxima- 
tion when studying non-axisymmetric models, be- 
cause rotating non-axisymmetric components (like 
bar or spirals) can intro d uce resonances (see e.g. 
Patsis fc Grosbol I (|l996f ). IPatsis et al.l (|2002h ). 



Our evolution calculations are presented to sim- 
ulate more realistic old NS distribution under the 
two-component Maxwellian initial random veloc- 
ity. We model the gravitational potential of the 
Galaxy following P90: 



sph 



disk 



halo 



(2) 



where $ sp h, 'I'haio and <i>disk are the spheroid, halo 
and disk components, respectively. 

For the spheroid and disk components one has: 



#i(R,z) = - 



GM; 



VR 2 + [a i + (z 2 + bj ! ) 1 / 2 ] 2 



(3) 



where R = \/x 2 + y 2 is the distance from the 
Galactic rotation axis and z is the distance from 
the Galactic disk plane. The subscript "i" repre- 
sents "sph" and "disk" . The values for the param- 
eters are taken from P90 and for the spheroid com- 
ponent they are: a sp h = 0.0 kpc, b sp h = 0.28 kpc 



2 



and M sph = 1.12 x lO 10 M Q ; while for the disk 
component: adisk = 3.7 kpc, bdisk = 0.20 kpc and 
M disk = 8.01 x 10 10 M . 

The halo component of the Galactic gravita- 
tional potential is: 



GM 



halo 



halo 



Vr 2 " 



\ In I 



arctan 



R 2 



VR2 



(4) 
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where r c = 6.0 kpc and M ha io = 5.0 x 10 lu M G 

2.3. NS initial position distribution 

It is generally accepted that the galactic 
z-distribution of mass ive objects is ap proxi- 



mately exponential ijBinnev fc Merrifield 1998: 



Mdzinarishvili fc Melikidze 1 12004) . This kind of 



the distribution can be theoretically explained by 
considering the dynamic equilibrium within the 
Galaxy. The initial z probability density func- 
tion of NSs in th e Galaxy has been proposed by 
iGott et al.l (1 1 9 7 Clf) and adopted by m any authors 
since then (e.g. iGonthier et al.l (|2002l )): 



Pz(z) 



1 

2hl 



exp 



where h z = 0.07 kpc is the scale height and: 



dz = 1 



,oo j 


-Nil 


/ TT exp 
Jo h z 





(5) 



(6) 



For the initial radial probability density func- 
tio n of the NSs we adopt t he same expression 



Arzoumanian et al.l (|2002h . As in P90, it tol- 



as 

lows a gamma function r(2,4.5), but has a radial 
outer boundary at 15 kpc rather than at 20 
This is motivated by the radial distribution of NS 
progenitors, i.e. population I massive stars. Al- 
though the Galaxy is believed to have a stellar 
disk ~ 15 kpc and a gaseous disk 15 ~ 25 kpc, 
NS progenitors hardly form in the gaseous disk, 
due to the considerable decreas e of the gas den- 
sity (jJones fc Lambourne 1120041) . 



1 Because of the rapid decrease of the Gamma func- 
tion with R we do not expect this modification to 
have a large impact on the results. 



The initial radial probability density function 
that we use is the following: 



R 

Pr(R) = a r-^— exp [-R/R 



cxpj 



where 



a r 



1 — exp 



1 



R 



exp 



(7) 



(8) 



We use i? cxp = 4.5 kpc, which gives a r ~ 1.183. 
The probability distribution is normalized to 1 
within the considered radial domain, i.e. from 
to 15 kpc. 

2.4. NS initial velocity distribution and 
equations of motion 

The NS initial velocity is calculated as the 
vector addition of three different velocities: (1) a 
Maxwcllian distribution, (2) a constant kick, and 
(3) the circular rotation velocities at the birth- 
place. 

Maxwellian distributions are usually used to 
represent the observed distribution of pulsar ve- 
locities. In this work we choose a two-component 
Maxwcllian distribution. One component includes 
40% of all NSs and has a velocity dispersion a v ~ 
90 km/s. The other one includes the remaining 
60% and has a v ~ 500 km/s, as proposed by 



Arzoumanian et al. I d2002h 



While for the kick velocity we adopt the con- 
ventional value of about 400 km/s for every sin- 
gle ob ject ( Hansen fc Phinnev~lll997 : Hobbs et al 
20051 ) . 



The initial circular rotation velocity of the NS 
is determined by 



d$\ 
rc=|R dR) 



1/2 



(9) 



where <& is the P90's gravitational potential in 
Eq. ©. 

The differential equations that describe the NS 
motion in the Galaxy can be expressed in the com- 
pact vector form as 



f = -V$ ( Vx' 



y 2 >z 



(10) 



where r= -\/x 2 + y 2 + z 2 is the spherical distance 
from the galactic center. NS orbits are numcr- 
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ically integrated with the fourth-order Runge- 
Kutta method. 

The NS total energy integral is used to control 
the accuracy of the integrations and in our simu- 
lations the total energy change is less than 1 part 
in 10 6 . The accuracy changes for different orbits, 
and generally simpler orbits are more accurate. 

3. NS orbits in the Galaxy 

The Poincare section technique is a way of 
presenting a trajectory in (n)-dimensional phase 
space in an (n— l)-dimcnsional space. By picking 
one phase element constant and plotting the values 
of the other elements each time the selected ele- 
ment has the desired value, an intersection surface 
is obtained. This technique has been used by sev- 
eral authors (e.g. Katsanikas fc Patsis ( 201 lh ) to 
analyze the structure of phase space in the neigh- 
borhood of stable periodic orbits in a 3D potential, 
and the properties of the invariant tori in the 4D 
spaces of section, under different galactic poten- 
tials. We use it here to study the 3-D NS trajec- 
tories and their 2-D projections. 

We plot the Poincare section for x > 0, and 
we fix y = to investigate the dynamical 3-D or- 
bits of NSs, as illustrated in row C in Fig. [TJ with 
varying the initial parameters. The phase space 
of NS's motion is 6-D, but since the total energy 
and angular momentum are conserved it is in fact 
only 4-D and its Poincare section is 3-D. The ini- 
tial condition (x,y, z,v x ,v y ,v z ) is reported under 
the corresponding column in the Figure. The NSs' 
motions are very diversified. 
(|2010bl) 



Wei et al 



investigated the gravita- 
tional potential of the Galactic disk and orbits 
of stars, and found that all of the orbits are sym- 
metric with respect to the galactic plane. Here 
we use the P90 gravitational potential and find 
that there are some non-symmetric orbits, see row 
D, columns /3 and 7 in Fig. Q] In row A we can 
see that when the motion range in the vertical 
direction becomes larger than the one in the ra- 
dial direction, the orbits become more irregular. 
See also same behavior for the projection on x-y 
plane in row B. While from row E, we see that 
the intersection points distribute in some regular 
lines on the projection of the Poincare section, 
which is essentially a closed curve. As such, the 



motion appears as a quasipcriodic orbit. However 
if the motion were exactly periodic, we would ex- 
pect that after some time, the star should return 
back to the same intersection point on the sur- 
face section, and this is not always the case in our 
simulations. 

According to P90, the dynamical behavior of 
NS populations is insensitive to the initial scale- 
height of progenitors. Here we see that NS orbits 
are like those of their progenitors: they are all 
basically rotating around the Galactic center, at 
different radial distance and uniformly. 

4. Simulation Results and Discussions 

In our calculation we obtain that the NS dis- 
tribution is steady after 10 9 yr. After this time 
we see that the NSs have greatly expanded in the 
radial direction, with the majority of them being 
located beyond 25 kpc from the Galactic center. 
More precisely 80% of the old NSs remains within 
25 kpc from the Galactic rotation axis (i.e. R < 25 
kpc), and 18% instead remains within 25 kpc from 
the Galactic center (i.e. r < 25 kpc). NSs moving 
in and out of the above range are in a dynamically 
equilibrium state. 

4.1. Radial distribution of NSs 

As we mentioned earlier, NSs are born in the 
region 0—15 kpc and later on they sp read to all 
radii. We follow an approach similar to IWei et al 



(|2010bf ) in order to investigate the characteristics 
of old NS distribution under the two-component 
Maxwellian initial random velocity. The normal- 
ized position probability density function that we 
find is shown in Fig. [5] We find that the distri- 
bution deviates from the initial distribution, i.e. 
r(2,4.5), due to the NS motion in the Galactic 
gravitational field. The distribution peak is now 
closer to the Galactic center. At first we fit the 
normalized R probability density function with 
the Gamma function T(a, A) as: 



R a— 

p R (R) = exp 



-R/A 



(11) 



The best fitting Gamma function is T(1.7, 5.2). 
The peak location of a generic Gamma function 
r(a, A) is at r p = A(a — 1). Using this expression 
for the initial radial distribution we get 4.5 kpc, 
while for the simulated distribution the peak is at 
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x=-B.3 90kpc v 1 , = 16.fl6km/s 
Y=-0.183kpc v y =-2Q1 .17km/s 
z = . HSkpc v = =-llG.81km/s 



x=-2 .234kpc v„=112. 65kra/s 
y=-2. 950kpc v y =-217 . B lkra/ s 
z=0.14 3kpc v E =131. 62kra/s 



x=3. 64 7kpc v x = 13 . B9kin/s 

y=-2 .319kpc v y = -67 . llkm/s 
z=-0.00 9kpc v E = -lB6.71km/: 



Fig. 1. — 3-D orbits of NSs. Panels in row A: trajectories; panels in row B: projections of the trajectories 
on the galactic plane; panels in row C: Poincare section x > 0, y = of the 3-D trajectory; panels in row D: 
projections on x-z plane of the Poincare sections; panels in row E: projections of the Poincare sections on 
x-v x plane. 



3.71 ± 0.17 kpc. The fitting results arc listed in 
Tables HJ Hand |3] 

Due to the unsatisfactory fitting of Gamma 
function, specially at the peak, we use the "ze ro 



point c orrected T(a, A)" proposed by IWei et al 
( 2010bf ). They modify the Gamma function 
r(a, A) by adding a constant A : 



p R (R) = A + A 



R" 



A' 



■ exp 



— R/A 



(12) 



If the points are close to a Gamma distribution 
function then the scatter will be small relative to 
the total variation in the values of the response 
variable. We adopt the coefficient of determina- 
tion (COD, also known as r-squared) to measure 
the fit quality. The closer COD to 1, the better the 
fit. Figure [2] indicates that the Gamma distribu- 
tion function is quite satisfactory with a COD of 
0.99, and Eq. (fl~2)) is acceptable for the case with 
the relative standard errors of the fitting parame- 
ters less than 5%. 

As for the T function, also for the "zero point 
corrected T" the peak position is at A(a — 1) and 
for the best fit case it is at 3.73 ± 0.19 kpc. We 
notice that the evolution of the NSs in the Galac- 
tic gravitational field makes the value of their R- 
distribution at R = kpc not exactly equal to 0. 
In other words, there is a "zero shift" , which is the 
total effect of the NS orbits shown in Sec. [3] 

4.2. Vertical distribution of NSs 

We consider the vertical distribution of the 
bound NSs in the whole Galactic disk with R < 
25 kpc and we find that it is not well described 
by a single exponential decay. For this reason we 
employ a double exponential profile: 

p z (z) = A x g(z) + A 1 x exp [-z/hi] 

+ A 2 X exp [-z/h 2 ] , (13) 

where A x g(z) represents the disk component, 
g(z) is step function which is 1 in the disk and 
outside, and hi and h 2 are the height scales of 
the two exponential contributions. Without loss 
of generality we can assume hi < h 2 and refer to 
Ai x exp [— z/hi] as the low-scale-height compo- 
nent and to A 2 x exp [— z/h 2 ] as the high-scale- 
hcight component. The probabilities for the low- 
scale-height and the high-scale-height component 



simulated data points 
standard Gamma function 
zero point corrected T(k, a) 




5 10 15 20 25 

Fig. 2. — Radial probability density distribu- 
tion. The crosses (black) indicate simulated data 
points; the dotted (green) line indicates the initial 
Gamma function; the two solid lines indicate the 
best fit with the "zero point corrected T(a, A)" 
(blue), and with the standard Gamma function 
(brown). Details about the fits are given in Ta- 
bles [TJ [2] and [3l respectively 
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are respectively: 

/>oo 

Pi— Ai X exp [-z/hi] dz = Ai X hi (14) 
Jo 

/>OC 

P 2 = / A 2 x exp [-z/h 2 ] dz = A 2 x h 2 . (15) 
Jo 

We also study the half density scale height of 
the disk zi/ 2 , defined as the height at which the to- 
tal probability density drops to 50% of the Galac- 
tic plane one. These results are shown in Table [T] 
We get a COD ~ 0.99 and relative standard errors 
< 1%, except that the relative standard error of 
p lies in the range 6.4 %. 



Table 1: Parameters of two exponential decay 
model Eq. ([TBI . 



parameter 


value 


relative error % 


A (kpc- 1 ) 


1.8 x 10" 5 


6.4 


Aitkpc- 1 ) 


1.87 


0.02 


hi (kpc) 


20.6 x 10~ 3 


0.04 


A 2 (kpc" 1 ) 


35.6 x 10~ 3 


0.09 


h 2 (kpc) 


1.55 


0.08 


COD 


0.999 




P2/P1 


1.45 


0.003 


zi/ 2 (kpc) 


17.6 x 10~ 3 


0.01 



Table 3: Best fit parameters for the "zero point 
corrected T(a, A)" . 



parameter 


value 


error % 


A (kpc- 1 ) 


5 x 10" 3 


6.3 


A 


95.6 x 10~ 3 


1.2 


a 


1.83 


0.58 


A (kpc) 


4.48 


1.1 


COD 


0.998 





Table 4: Best fit parameters for linear relations of 
the radial distributions of the two scale heights hi 
and h 2 .^ is for fitting the region inward of 4.25 kpc 
and * is for fitting the region outward. 



parameter 


value 


error % 


ki 


13 x 10~ 3 


0.26 


bi 


12.8 x 10~ 3 


0.005 


k 2 


18.4 x 10~ 3 


0.0007 


b 2 


0.03 


0.002 


k* 


0.05 


0.66 


b* 


0.65 


0.16 



Table 2: Best fit results for the standard Gamma 
function. 



parameter 


value 


relative error % 


A 


1.13 


0.49 


a 


1.71 


0.54 


A(kpc) 


5.21 


0.86 


COD 


0.996 





Observational studies of the Galactic disk re- 
veal that it can be well described in terms of 
the two components model with a t hin disk and 
a thick disk component (se e e.g. 



Chen et al. 



(|200lh . lKaempf et al. ll|2005h . lWei et al.1 (l2010a|) 1 



Our simulation confirms the validity of this model 



and shows that the z hierarchy effect can be re- 
garded as the result of the dynamical evolution of 
the old NSs originated from the Galactic disk. 

4.3. Scale-Height vs R relation 

We now consider the z distribution at differ- 
ent Galactic radial distances R from the Galactic 
center. To this end, we divide R with 0.5 kpc 
spacing from kpc to 25 kpc and get 50 parts, 
then we analyze in details the NS z distribution in 
each part. The two exponential decay of Eq. ([13"j) 
is still employed to study the case in each spac- 
ing segment of R, to derive the high-scale-heights 
h 2 , low-scale-heights hi, and the ratios of the two 
components P 2 /Pi. 

The results are shown in Fig. |3] In each seg- 
ment of R the two exponential model is still sig- 
nificantly effective. As a result, the CODs of both 
distributions in Fig. [3] is quite similar ~ 0.99, and 
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with very small relative standard errors < 0.01 . 
The fitting results are listed in Table |U 

The relationship between the two scale heights 
and R can be described by the following linear 
model. 

h(R) = kxR + b, (16) 

We also plot the radial dependence of the half den- 
sity scale heights in Fig. [3] The fitting results of 
parameters k and b are listed in Table SJ 

The low-scale-heights can be depicted by a lin- 
ear function of R within the whole range of the 
Galactic disk — 25 kpc. The slope of the fit- 
ting line is small, which means that the changing 
of this component within — 25 kpc of R is not 
large. 

For the high-scale-height components, there ex- 
ists a point at R = 4.25 kpc, where the behavior 
changes. Both sides of this point have a linear ra- 
dial dependence but with different slopes: the one 
inside R < 4.25 kpc is smaller than the one out- 
side R > 4.25 kpc. We notice that R = 4.25 kpc 
corresponds to the observed Rq of the HI disk. 

As regards the R distribution of P 2 /Pi, we 
can see the difference o f old N S distribution un- 
der lArzoumanian et al. I (120021) from those under 
iHobbs et al.l (120051) andlFaucher-Giguere fc Kaspi 
(|2006l ) in I Wei et al.l (j2010b|). Where the R distri- 
butions of P2/P1 have three distinct parts clearly, 
and the high-scale-height distributions have not 
points like the observed Ro of the HI disk. The 
ratio of the two components leads to the increase 
of growth slowly and smoothly with R in the whole 
Galaxy. 

Table 5: Best fit parameters for linear relations of 
the radial distributions of the half density of two 
scale heights Eq. (|T6|) . 



parameter 


value 


error % 


k 


14.9 x 10~ 3 


0.01 


1) 


10.5 x 10- 3 


0.01 




Fig. 3. — Scale-heights in the double exponential 
decay model see Eq. (fl"3"|) . Top panel: ratio of 
P2/P1. Bottom panel: radial distributions of the 
two scale heights hi and h.2 . For both of the panels 
the crosses indicate simulated data points, while 
the solid lines indicate the best fit line. The fitting 
parameters are given in Table 21 
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Fig. 4. — Radial distributions of the half density 
scale heights in the region R< 25 kpc with the 
fitting parameters of Eq. (|16|) for the straight lines. 
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The NS scale-heights generally increase from 
the Galactic center to the edge of Galactic disk. 
This phenomenon is independent of the initial ve- 
locity distribution of the NSs, and it is due to the 
action of the Galactic gravitational field on the 
NSs. Our calculation shows that the heights of 
the orbits of the NSs generally decrease towards 
the Galactic center (see Sec. [3]). P90 calculated 
the half density scale height in the vicinity of the 
Sun (R = 8 kpc) and obtained zi/ 2 = 0.2 kpc. 
In our present updated version the corresponding 
value is 0.2 ± 0.0004 kpc. 

5. Discussions and Conclusions 

In this paper we have investigated the space 
and velocity distribution of old neutron stars 
(NSs) in our Galaxy. We assume that the ini- 
tial velocity distribution is the result of three 
components: the kick velocity, the circular ve- 
locity and the M a xwelli an velocity (following 
" (2Q0l). 



Arzoumanian et al 



For the initial posi- 
tion distribution instead we assume that it follows 
a r function. As regards the Galactic gravita - 



tional potential, we follow the iPaczvhski 



prescription, which is of course only and approx- 
imation, since it does not take into account any 
inhomogeneity within the Galaxy. However this 
is suitable for the simplified analysis that we are 
developing here. We have then used Monte Carlo 
simulations to let the NSs evolve and have shown 
3-D NS orbits and Poincare sections of the phase 
space. 

It is evident that the irregular character of the 
motion of NSs increases when the vertical direc- 
tion becomes larger than radial direction. Another 
remarkable finding is that there are some signifi- 
cant diffractions in the symmetric of the orbits, 
which may effects of supcrnovae kicks. 

Our numerical results show that NSs have a 
very large radial Galactic expansion. The major- 
ity of them (80%) falls within 25 kpc from the 
Galactic rotation axis (R < 25 kpc), and 18% 
instead remains within 25 kpc from the Galactic 
center (r < 25 kpc). An important aspect is that 
the total number of NSs moving in and out of the 
above range is in a dynamically equilibrium state 
after 10 9 yr. 

The radial probability density distribution de- 
viates from the initial distribution, and has a peak 



which is closer to the Galactic center. The analy- 
sis of the vertical and radial distributions clearly 
show that the orbits of the NSs decrease toward 
the Galactic center within different scale heights. 

Qualitatively, the simulated old NSs disk, es- 
pecially the middle and outer components, keeps 
the observed HI disk in moderation. Although the 
old NSs and their progenitors have different radial 
and vertical distribution, we find that the shapes 
of their orbits are quite similar in the HI clouds 
regions. 

The results of this work will constitute the base 
for further studies on NS properties. Such re- 
search could be helpful for the detection of old 
NSs via their gravitational microlensing that re- 
sult in the variation in the brig htness of t h e dis- 
tant active galactic nuclei (e.g. Hawkins ( 19931 
2002])). Another way for detecting the old NSs 



is throu gh the interaction with the interstellar 
medium (IPopov et al. I l2000h . 



As subsequent steps we plan to (1) apply the 
three exponential decay model in studying the NS 
vertical distribution with more detail and (2) to 
use different models of the Galactic potential to 
investigate specific parts of our Galaxy. 

Acknowledgments 

This work is supported by the National Natu- 
ral Science Foundation of China (NSFC 10773017, 
NSFC 10773034) and National Basic Research 
Program of China (2009CB824800, 2012CB821800). 
Chinese Academy of Sciences and National Astro- 
nomical Observatory of China (NAOC) of CAS 
has supported this work by the Silk Road Project 
(CAS Grant Number 2009S1-5). L.N. is currently 
supported by a Chinese Academy of Sciences fel- 
lowship for young international scientists (Grant 
Number 2010Y2JB12). 

REFERENCES 

Arzoumanian, Z., Chernoff, D. F., & Cordes, J. M. 
Astrophys. J. 568, 289 (2002) 

Bhattacharya D. & van den Heuvel E. P. J., Phys. 
Rep., 203, 1 (1991) 

Binney J. & Merrifield M., Galactic Astronomy 
(Princeton: Princeton University Press) (1998) 



10 



Caldwell J. & Ostriker J., Astrophys. J. 251, 61 
(1981) 

Chen B., Stoughton C, Smith J. A., ct al. Astro- 
phys. J. 553, 184 (2001) 

Cordcs J. M., & ChcrnofTD. F., Astrophys. J. 505, 
315 (1998) 



Cordes, 



M., & Lazio, 



T. 



W. 



arXiv:astro-ph/0207156 (2002) 



Faucher-Giguere C.-A. & Kaspi V. M., Astrophys. 
J. 643, 332 (2006) 

Frei Z., Huang X. & Paczyhski B., ApJ. 643, 332 
(1992) 

Gonthier P. L., Ouellette M. S., Berrier J. et al. 
Astrophys. J. 565, 482 (2002) 

Gott J. R., Gunn J. E., & Ostriker J. P., Astro- 
phys. J. 160, L91 (1970) 

Hansen B. M. S. & Phinncy E. S., Mon. Not. R. 
Astron. Soc. 291, 569 (1997) 

Hawkins M.R.S., Nature 366, 242 (1993) 

Hawkins M.R.S., Mon. Not. R. Astron. Soc. 329, 
76 (2002) 

Hobbs G., Lorimer D. R., Lyne A. G., & Kramer 
M., Mon. Not. R.Astron. Soc. 360, 974 (2005) 

Jones M. H, & Lambourne R.J. A. An Introduction 
to Galaxies and Cosmology. Cambridge: Cam- 
bridge University Press, 7 (2004) 

Kaempf T. A., de Boer K. S., & Altmann M., As- 
tron. & Astrophys, 432, 879 (2005) 

Katsanikas M. & Patsis P. A., Int. J.Bif. Chaos, 
21, 467 (2011) 

Lorimer D. R. Living Rev. Relativity, 11,8 (2008) 



Mdzinarishvili T. G. & Melikidze G. I., Astron.fe 
Astrophys., 425, 1009 (2004) 

Ofek E. O., PASP, 121, 814 (2009) 

Paczyhski B., Astrophys. J. 348, 485 (1990) 

Patsis P. A. & Grosb0l P., Astron.fe Astrophys., 
315, 371 (1996) 

Patsis P. A., Athanassoula E., Grosb0l P. et al. 
Mon. Not. R. Astron. Soc. 355, 1049 (2002) 

Perets H. B, Wu X., Zhao H. S. et al. ApJ, 697, 
2097 (2009) 

Popov S. B., Colpi M., Treves A. et al. ApJ, 530, 
896 (2000) 

Sartore N., Ripamonti E., Treves A. & Turolla R., 
Astron.fe Astrophys., 510, A23 (2010) 

Taylor J. H., & Cordes J. M., Astrophys. J. 411, 
674 (1993) 

Wei Y. C, Taani A., Pan Y. Y. et al., Chin. Phys. 
Lett., 27, 9801 (2010a) 

Wei Y. C, Chengmin C. M., Xinji W. et al., Scince 
in China, 53, 1939 (2010b) 



http://relativity.livingreviews.org/Articles/lrr-2008-8/ 



Lyne A. G., Anderson B., & Salter M. J., Mon. 
Not. R. Astron. Soc. 201, 503 (1982) 

Lyne A. G., Burgay M., Kramer M. et al. Science, 
303, 1153 (2004) 

Lyne A. G. & Graham-Smith, F., Pulsar As- 
tronomy, Cambridge Astrophysics Series, Cam- 
bridge University Press, 3ed Edit. (2007) 



This 2-column preprint was prepared with the AAS IATjrjX 
macros v5.2. 



11 



